Go Big or Go Home Part 2 - Working and Visualising on cluster

In this notebook we finally do our larger area. We're going to need some better visualisation tools and it would be great not to bring the results back to the Jupyter notebook but to leverage the dask clusters resources during visualisation. We'll be using some dask-aware visualiation libraries (holoviews and datashader) to do the heavy lifting.

Let's begin by starting up our cluster and sizing it appropriately to our computational task.

Time to go big!¶

All the code here is the same as the conclusion from the previous notebook, except we'll make the cluster bigger with 10 workers instead of 4. We'll also make the masking and NDVI calculation into a python function since we won't be making any changes to that now.

We'll use the same ROI and time period for this run and we're using all the techniques so far to reduce the computation time:

  1. Dask chunk size selection
  2. Only loading the measurements we intend on using in this calculation to save on the task graph optimisation time
In [27]:
# Initialize the Gateway client
from dask.distributed import Client
from dask_gateway import Gateway

number_of_workers = 10 

gateway = Gateway()

clusters = gateway.list_clusters()
if not clusters:
    print('Creating new cluster. Please wait for this to finish.')
    cluster = gateway.new_cluster()
else:
    print(f'An existing cluster was found. Connecting to: {clusters[0].name}')
    cluster=gateway.connect(clusters[0].name)

cluster.scale(number_of_workers)

client = cluster.get_client()
client

An existing cluster was found. Connecting to: easihub.b444f7be6e5f4b9888259e752ac0d2b1
Out[27]:

Client

Client-ad70b3b6-064b-11ee-82af-8e60251f35b2

Connection method: Cluster object Cluster type: dask_gateway.GatewayCluster
Dashboard: https://hub.adias.aquawatchaus.space/services/dask-gateway/clusters/easihub.b444f7be6e5f4b9888259e752ac0d2b1/status

Cluster Info

GatewayCluster

  • Name: easihub.b444f7be6e5f4b9888259e752ac0d2b1
  • Dashboard: https://hub.adias.aquawatchaus.space/services/dask-gateway/clusters/easihub.b444f7be6e5f4b9888259e752ac0d2b1/status
In [2]:
import pyproj
pyproj.set_use_global_context(True)

import git
import sys, os
from dateutil.parser import parse
from dateutil.relativedelta import relativedelta
from dask.distributed import Client, LocalCluster
import datacube
from datacube.utils import masking
from datacube.utils.aws import configure_s3_access

# EASI defaults
os.environ['USE_PYGEOS'] = '0'
repo = git.Repo('.', search_parent_directories=True).working_tree_dir
if repo not in sys.path: sys.path.append(repo)
from easi_tools import EasiDefaults, notebook_utils
easi = EasiDefaults()

Successfully found configuration for deployment "adias"
In [3]:
dc = datacube.Datacube()
configure_s3_access(aws_unsigned=False, requester_pays=True, client=client)
Out[3]:
<botocore.credentials.Credentials at 0x7f6d148e0940>
In [4]:
# Get the centroid of the coordinates of the default extents
central_lat = sum(easi.latitude)/2
central_lon = sum(easi.longitude)/2
# central_lat = -42.019
# central_lon = 146.615

# Set the buffer to load around the central coordinates
# This is a radial distance for the bbox to actual area so bbox 2x buffer in both dimensions
buffer = 0.8

# Compute the bounding box for the study area
study_area_lat = (central_lat - buffer, central_lat + buffer)
study_area_lon = (central_lon - buffer, central_lon + buffer)

# Data product
products = easi.product('landsat')

# Set the date range to load data over
set_time = easi.time
set_time = (set_time[0], parse(set_time[0]) + relativedelta(years=1))
# set_time = ("2021-01-01", "2021-12-31")

# Selected measurement names (used in this notebook)
alias = easi.aliases('landsat')
measurements = [alias[x] for x in ['qa_band', 'red', 'nir']]

# Set the QA band name and mask values
qa_band = alias['qa_band']
qa_mask = easi.qa_mask('landsat')

# Set the resampling method for the bands
resampling = {qa_band: "nearest", "*": "average"}

# Set the coordinate reference system and output resolution
set_crs = easi.crs('landsat')  # If defined, else None
set_resolution = easi.resolution('landsat')  # If defined, else None
# set_crs = "epsg:3577"
# set_resolution = (-30, 30)

# Set the scene group_by method
group_by = "solar_day"
In [5]:
def masked_seasonal_ndvi(dataset):
    # Identify pixels that are either "valid", "water" or "snow"
    cloud_free_mask = masking.make_mask(dataset[qa_band], **qa_mask)
    # Apply the mask
    cloud_free = dataset.where(cloud_free_mask)

    # Calculate the components that make up the NDVI calculation
    band_diff = cloud_free[alias['nir']] - cloud_free[alias['red']]
    band_sum = cloud_free[alias['nir']] + cloud_free[alias['red']]
    # Calculate NDVI
    ndvi = None
    ndvi = band_diff / band_sum

    return ndvi.groupby("time.season").mean("time")  # Calculate the seasonal mean

dataset = None # clear results from any previous runs
dataset = dc.load(
            product=products,
            x=study_area_lon,
            y=study_area_lat,
            time=set_time,
            measurements=measurements,
            resampling=resampling,
            output_crs=set_crs,
            resolution=set_resolution,
            dask_chunks = {"time":2, "x":3072, "y":3072},
            group_by=group_by,
        )

ndvi_unweighted = masked_seasonal_ndvi(dataset)
In [6]:
print(f"dataset size (GiB) {dataset.nbytes / 2**30:.2f}")
print(f"ndvi_unweighted size (GiB) {ndvi_unweighted.nbytes / 2**30:.2f}")

dataset size (GiB) 37.06
ndvi_unweighted size (GiB) 2.86
In [7]:
client.wait_for_workers(n_workers=10)
In [8]:
cluster

VBox(children=(HTML(value='<h2>GatewayCluster</h2>'), HBox(children=(HTML(value='\n<div>\n<style scoped>\n    …
In [9]:
%%time
actual_result = ndvi_unweighted.compute()

CPU times: user 4.34 s, sys: 1.88 s, total: 6.22 s
Wall time: 3min 50s

You'll notice the computation time is slightly faster with more workers - we're IO bound so more workers means more available IO bandwidth and threads. It's not 2-3 x faster though - we're wasting a lot of resources because we can't actually use all of that extra power.

Tip: More isn't always better. Be mindful of your computational resource usage and cost. This size cluster is a tremendous waste for this size computational job. Size things appropriately.

And visualise the result

In [10]:
actual_result.sel(season='DJF').plot()
Out[10]:
<matplotlib.collections.QuadMesh at 0x7f6d143c78e0>

We'll save the coordinates of this section from the array (as slices) so we can use them later for visualising the same ROI from a larger dataset.

In [11]:
x_slice = slice(ndvi_unweighted.x[0], ndvi_unweighted.x[-1])
y_slice = slice(ndvi_unweighted.y[0], ndvi_unweighted.y[-1])

Now for a bigger area¶

Let's change the area extent to about 4 degrees square.

In [12]:
# Compute the bounding box for the study area
buffer = 2
# Compute the bounding box for the study area
study_area_lat = (central_lat - buffer, central_lat + buffer)
study_area_lon = (central_lon - buffer, central_lon + buffer)
In [13]:
# Check the map below to see if you are including the area that you want. For this example, it would be best to not include too much water.
from dea_tools.plotting import display_map
display_map(study_area_lon, study_area_lat)

/env/lib/python3.10/site-packages/dea_tools/plotting.py:313: FutureWarning: This function is deprecated. See: https://pyproj4.github.io/pyproj/stable/gotchas.html#upgrading-to-pyproj-2-from-pyproj-1
  all_longitude, all_latitude = transform(Proj(crs), Proj('EPSG:4326'), all_x,
Out[13]:
Make this Notebook Trusted to load map: File -> Trust Notebook
In [14]:
dataset = None # clear results from any previous runs
dataset = dc.load(
            product=products,
            x=study_area_lon,
            y=study_area_lat,
            time=set_time,
            measurements=measurements,
            resampling=resampling,
            output_crs=set_crs,
            resolution=set_resolution,
            dask_chunks = {"time":2, "x":3072, "y":3072},
            group_by=group_by,
        )

ndvi_unweighted = masked_seasonal_ndvi(dataset)

Before we compute anything let's take a look at our result's shape and size

In [15]:
print(f"dataset size (GiB) {dataset.nbytes / 2**30:.2f}")
print(f"ndvi_unweighted size (GiB) {ndvi_unweighted.nbytes / 2**30:.2f}")

dataset size (GiB) 375.76
ndvi_unweighted size (GiB) 17.89

This is now much bigger!

The result is getting on the large size for the notebook node so we will need to pay attention to data locality and the size of results being processed. The cluster has a LOT more memory than the notebook node; bring too much back to the notebook and the notebook will crash.

Tip: Be mindful of the size of the results and their data locality.

Now let's check the shape, tasks and chunks

In [16]:
dataset
Out[16]:
<xarray.Dataset>
Dimensions:      (time: 112, y: 24512, x: 24494)
Coordinates:
  * time         (time) datetime64[ns] 2022-02-04T18:45:01.040642 ... 2023-01...
  * y            (y) float64 1.413e+07 1.413e+07 ... 1.339e+07 1.339e+07
  * x            (x) float64 5.1e+06 5.1e+06 5.1e+06 ... 5.835e+06 5.835e+06
    spatial_ref  int32 32650
Data variables:
    qa_pixel     (time, y, x) uint16 dask.array<chunksize=(2, 3072, 3072), meta=np.ndarray>
    red          (time, y, x) uint16 dask.array<chunksize=(2, 3072, 3072), meta=np.ndarray>
    nir08        (time, y, x) uint16 dask.array<chunksize=(2, 3072, 3072), meta=np.ndarray>
Attributes:
    crs:           epsg:32650
    grid_mapping:  spatial_ref
xarray.Dataset
    • time: 112
    • y: 24512
    • x: 24494
    • time
      (time)
      datetime64[ns]
      2022-02-04T18:45:01.040642 ... 2...
      units :
      seconds since 1970-01-01 00:00:00
      array(['2022-02-04T18:45:01.040642000', '2022-02-06T18:32:39.001148000',
       '2022-02-11T18:51:09.926245000', '2022-02-13T18:38:47.450521000',
       '2022-02-15T18:27:12.621325000', '2022-02-20T18:44:54.824298000',
       '2022-02-22T18:32:33.208594000', '2022-02-27T18:51:05.324304000',
       '2022-03-01T18:38:43.403434000', '2022-03-03T18:27:09.162918000',
       '2022-03-08T18:44:52.352548000', '2022-03-10T18:32:29.902138000',
       '2022-03-15T18:50:59.645542000', '2022-03-17T18:38:36.674111000',
       '2022-03-19T18:27:01.370411000', '2022-03-24T18:44:41.976128000',
       '2022-03-26T18:32:18.534534000', '2022-03-31T18:50:48.087361000',
       '2022-04-02T18:38:26.196010000', '2022-04-04T18:26:53.904150000',
       '2022-04-09T18:44:41.660779000', '2022-04-11T18:32:20.904923000',
       '2022-04-16T18:50:53.974035000', '2022-04-18T18:38:31.913586000',
       '2022-04-20T18:26:57.478599000', '2022-04-25T18:44:39.910946000',
       '2022-04-27T18:32:17.063038000', '2022-05-02T18:50:52.394776000',
       '2022-05-04T18:38:32.304289000', '2022-05-06T18:26:59.822251000',
       '2022-05-11T18:44:47.359594000', '2022-05-13T18:32:26.621212000',
       '2022-05-18T18:51:00.436500000', '2022-05-20T18:38:38.910059000',
       '2022-05-22T18:27:04.890363000', '2022-05-27T18:44:48.250613000',
       '2022-05-29T18:32:28.727749000', '2022-06-03T18:51:05.898078000',
       '2022-06-05T18:38:45.960055000', '2022-06-07T18:27:13.702518000',
       '2022-06-12T18:45:01.955709000', '2022-06-14T18:32:41.535010000',
       '2022-06-19T18:51:16.329680000', '2022-06-21T18:38:55.358354000',
       '2022-06-23T18:27:22.024943000', '2022-06-28T18:45:07.537274000',
       '2022-06-30T18:32:46.074409000', '2022-07-05T18:51:18.547117000',
       '2022-07-07T18:38:56.758064000', '2022-07-09T18:27:22.619818000',
       '2022-07-14T18:45:06.061210000', '2022-07-16T18:32:45.301782000',
       '2022-07-21T18:51:22.298885000', '2022-07-23T18:39:02.212774000',
       '2022-07-25T18:27:29.794387000', '2022-07-30T18:45:17.725644000',
       '2022-08-01T18:32:57.253868000', '2022-08-06T18:51:32.302026000',
       '2022-08-08T18:39:11.550933000', '2022-08-10T18:27:38.463263000',
       '2022-08-15T18:45:24.556392000', '2022-08-17T18:33:03.284411000',
       '2022-08-22T18:51:36.142715000', '2022-08-24T18:39:14.519012000',
       '2022-08-26T18:27:40.622198000', '2022-08-31T18:45:24.812984000',
       '2022-09-02T18:33:02.829160000', '2022-09-07T18:51:38.534099000',
       '2022-09-09T18:39:17.982204000', '2022-09-16T18:45:31.348066000',
       '2022-09-18T18:33:10.105949000', '2022-09-23T18:51:42.960724000',
       '2022-09-25T18:39:21.265191000', '2022-09-27T18:27:47.207724000',
       '2022-10-02T18:45:30.599018000', '2022-10-04T18:33:08.094669000',
       '2022-10-09T18:51:37.077805000', '2022-10-11T18:39:13.505811000',
       '2022-10-13T18:27:37.669230000', '2022-10-18T18:45:26.431834000',
       '2022-10-20T18:33:06.096236000', '2022-10-25T18:51:41.069181000',
       '2022-10-27T18:39:20.161786000', '2022-10-29T18:27:46.833208000',
       '2022-11-03T18:45:31.741431000', '2022-11-05T18:33:09.719987000',
       '2022-11-10T18:51:39.986856000', '2022-11-12T18:39:17.046216000',
       '2022-11-14T18:27:41.667336000', '2022-11-19T18:45:24.533861000',
       '2022-11-21T18:33:03.664615000', '2022-11-26T18:51:37.421591000',
       '2022-11-28T18:39:16.031024000', '2022-11-30T18:27:42.274155000',
       '2022-12-05T18:45:50.377394000', '2022-12-07T18:33:04.379724000',
       '2022-12-12T18:51:34.594981000', '2022-12-14T18:39:11.696450000',
       '2022-12-16T18:27:36.338531000', '2022-12-21T18:45:40.193300000',
       '2022-12-23T18:32:54.280347000', '2022-12-28T18:51:28.515857000',
       '2022-12-30T18:39:31.118057000', '2023-01-01T18:27:33.470203000',
       '2023-01-06T18:46:05.172439000', '2023-01-08T18:32:54.991488000',
       '2023-01-15T18:39:00.610083000', '2023-01-17T18:27:48.368126000',
       '2023-01-22T18:45:13.659001000', '2023-01-24T18:32:53.086260000',
       '2023-01-29T18:51:26.822968000', '2023-01-31T18:39:05.296909000'],
      dtype='datetime64[ns]')
    • y
      (y)
      float64
      1.413e+07 1.413e+07 ... 1.339e+07
      units :
      metre
      resolution :
      -30.0
      crs :
      epsg:32650
      array([14126265., 14126235., 14126205., ..., 13390995., 13390965., 13390935.])
    • x
      (x)
      float64
      5.1e+06 5.1e+06 ... 5.835e+06
      units :
      metre
      resolution :
      30.0
      crs :
      epsg:32650
      array([5099805., 5099835., 5099865., ..., 5834535., 5834565., 5834595.])
    • spatial_ref
      ()
      int32
      32650
      spatial_ref :
      PROJCS["WGS 84 / UTM zone 50N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",117],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32650"]]
      grid_mapping_name :
      transverse_mercator
      array(32650, dtype=int32)
    • qa_pixel
      (time, y, x)
      uint16
      dask.array<chunksize=(2, 3072, 3072), meta=np.ndarray>
      units :
      bit_index
      nodata :
      1
      flags_definition :
      {'snow': {'bits': 5, 'values': {'0': 'not_high_confidence', '1': 'high_confidence'}}, 'clear': {'bits': 6, 'values': {'0': 'not_clear', '1': 'clear'}}, 'cloud': {'bits': 3, 'values': {'0': 'not_high_confidence', '1': 'high_confidence'}}, 'water': {'bits': 7, 'values': {'0': 'land_or_cloud', '1': 'water'}}, 'cirrus': {'bits': 2, 'values': {'0': 'not_high_confidence', '1': 'high_confidence'}}, 'nodata': {'bits': 0, 'values': {'0': False, '1': True}}, 'qa_pixel': {'bits': [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15], 'values': {'1': 'Fill', '2': 'Dilated Cloud', '4': 'Cirrus', '8': 'Cloud', '16': 'Cloud Shadow', '32': 'Snow', '64': 'Clear', '128': 'Water', '256': 'Cloud Confidence low bit', '512': 'Cloud Confidence high bit', '1024': 'Cloud Shadow Confidence low bit', '2048': 'Cloud Shadow Confidence high bit', '4096': 'Snow Ice Confidence low bit', '8192': 'Snow Ice Confidence high bit', '16384': 'Cirrus Confidence low bit', '32768': 'Cirrus Confidence high bit'}, 'description': 'Level 2 pixel quality'}, 'cloud_shadow': {'bits': 4, 'values': {'0': 'not_high_confidence', '1': 'high_confidence'}}, 'dilated_cloud': {'bits': 1, 'values': {'0': 'not_dilated', '1': 'dilated'}}, 'cloud_confidence': {'bits': [8, 9], 'values': {'0': 'none', '1': 'low', '2': 'medium', '3': 'high'}}, 'cirrus_confidence': {'bits': [14, 15], 'values': {'0': 'none', '1': 'low', '2': 'reserved', '3': 'high'}}, 'snow_ice_confidence': {'bits': [12, 13], 'values': {'0': 'none', '1': 'low', '2': 'reserved', '3': 'high'}}, 'cloud_shadow_confidence': {'bits': [10, 11], 'values': {'0': 'none', '1': 'low', '2': 'reserved', '3': 'high'}}}
      crs :
      epsg:32650
      grid_mapping :
      spatial_ref
      Array Chunk
      Bytes 125.25 GiB 36.00 MiB
      Shape (112, 24512, 24494) (2, 3072, 3072)
      Dask graph 3584 chunks in 2 graph layers
      Data type uint16 numpy.ndarray
      		24494 24512 112
    • red
      (time, y, x)
      uint16
      dask.array<chunksize=(2, 3072, 3072), meta=np.ndarray>
      units :
      reflectance
      nodata :
      0
      crs :
      epsg:32650
      grid_mapping :
      spatial_ref
      Array Chunk
      Bytes 125.25 GiB 36.00 MiB
      Shape (112, 24512, 24494) (2, 3072, 3072)
      Dask graph 3584 chunks in 2 graph layers
      Data type uint16 numpy.ndarray
      		24494 24512 112
    • nir08
      (time, y, x)
      uint16
      dask.array<chunksize=(2, 3072, 3072), meta=np.ndarray>
      units :
      reflectance
      nodata :
      0
      crs :
      epsg:32650
      grid_mapping :
      spatial_ref
      Array Chunk
      Bytes 125.25 GiB 36.00 MiB
      Shape (112, 24512, 24494) (2, 3072, 3072)
      Dask graph 3584 chunks in 2 graph layers
      Data type uint16 numpy.ndarray
      		24494 24512 112
    • time
      PandasIndex
      PandasIndex(DatetimeIndex(['2022-02-04 18:45:01.040642', '2022-02-06 18:32:39.001148',
               '2022-02-11 18:51:09.926245', '2022-02-13 18:38:47.450521',
               '2022-02-15 18:27:12.621325', '2022-02-20 18:44:54.824298',
               '2022-02-22 18:32:33.208594', '2022-02-27 18:51:05.324304',
               '2022-03-01 18:38:43.403434', '2022-03-03 18:27:09.162918',
               ...
               '2022-12-30 18:39:31.118057', '2023-01-01 18:27:33.470203',
               '2023-01-06 18:46:05.172439', '2023-01-08 18:32:54.991488',
               '2023-01-15 18:39:00.610083', '2023-01-17 18:27:48.368126',
               '2023-01-22 18:45:13.659001', '2023-01-24 18:32:53.086260',
               '2023-01-29 18:51:26.822968', '2023-01-31 18:39:05.296909'],
              dtype='datetime64[ns]', name='time', length=112, freq=None))
    • y
      PandasIndex
      PandasIndex(Float64Index([14126265.0, 14126235.0, 14126205.0, 14126175.0, 14126145.0,
              14126115.0, 14126085.0, 14126055.0, 14126025.0, 14125995.0,
              ...
              13391205.0, 13391175.0, 13391145.0, 13391115.0, 13391085.0,
              13391055.0, 13391025.0, 13390995.0, 13390965.0, 13390935.0],
             dtype='float64', name='y', length=24512))
    • x
      PandasIndex
      PandasIndex(Float64Index([5099805.0, 5099835.0, 5099865.0, 5099895.0, 5099925.0, 5099955.0,
              5099985.0, 5100015.0, 5100045.0, 5100075.0,
              ...
              5834325.0, 5834355.0, 5834385.0, 5834415.0, 5834445.0, 5834475.0,
              5834505.0, 5834535.0, 5834565.0, 5834595.0],
             dtype='float64', name='x', length=24494))
  • crs :
    epsg:32650
    grid_mapping :
    spatial_ref

Looking at the red data variable we can see about 50 GiB for the array, 36 MiB per chunk and 5526 tasks. Noting the nir and qa_band will be similarly shaped and size.

The number of tasks is climbing so we can expect an increase in task graph optimisation time.

Chunk size and tasks seems okay, but we will monitor the dask dashboard in case there are issues with temporaries causing workers to spill to disk if memory is too full.

The chunking is resulting in some slivers, particularly on the y axis. Let's modify the y chunk size so these slivers don't exist as its blowing out the tasks and is likely unnecessary. We will calculate the x and y chunk sizes below to get a nice fit. Make sure to check the chunk size afterwards to make sure it doesn't get too large. If it does we can make the chunks smaller to reduce slivers too.

In [17]:
from math import ceil

y_chunks = ceil(dataset.dims['y']/5)
y_chunks
Out[17]:
4903
In [18]:
dataset = None # clear results from any previous runs
dataset = dc.load(
            product=products,
            x=study_area_lon,
            y=study_area_lat,
            time=set_time,
            measurements=measurements,
            resampling=resampling,
            output_crs=set_crs,
            resolution=set_resolution,
            dask_chunks = {"time":2, "x":3072, "y":y_chunks},
            group_by=group_by,
        )

ndvi_unweighted = masked_seasonal_ndvi(dataset)

Now recheck our chunk size and tasks for red

In [19]:
dataset
Out[19]:
<xarray.Dataset>
Dimensions:      (time: 112, y: 24512, x: 24494)
Coordinates:
  * time         (time) datetime64[ns] 2022-02-04T18:45:01.040642 ... 2023-01...
  * y            (y) float64 1.413e+07 1.413e+07 ... 1.339e+07 1.339e+07
  * x            (x) float64 5.1e+06 5.1e+06 5.1e+06 ... 5.835e+06 5.835e+06
    spatial_ref  int32 32650
Data variables:
    qa_pixel     (time, y, x) uint16 dask.array<chunksize=(2, 4903, 3072), meta=np.ndarray>
    red          (time, y, x) uint16 dask.array<chunksize=(2, 4903, 3072), meta=np.ndarray>
    nir08        (time, y, x) uint16 dask.array<chunksize=(2, 4903, 3072), meta=np.ndarray>
Attributes:
    crs:           epsg:32650
    grid_mapping:  spatial_ref
xarray.Dataset
    • time: 112
    • y: 24512
    • x: 24494
    • time
      (time)
      datetime64[ns]
      2022-02-04T18:45:01.040642 ... 2...
      units :
      seconds since 1970-01-01 00:00:00
      array(['2022-02-04T18:45:01.040642000', '2022-02-06T18:32:39.001148000',
       '2022-02-11T18:51:09.926245000', '2022-02-13T18:38:47.450521000',
       '2022-02-15T18:27:12.621325000', '2022-02-20T18:44:54.824298000',
       '2022-02-22T18:32:33.208594000', '2022-02-27T18:51:05.324304000',
       '2022-03-01T18:38:43.403434000', '2022-03-03T18:27:09.162918000',
       '2022-03-08T18:44:52.352548000', '2022-03-10T18:32:29.902138000',
       '2022-03-15T18:50:59.645542000', '2022-03-17T18:38:36.674111000',
       '2022-03-19T18:27:01.370411000', '2022-03-24T18:44:41.976128000',
       '2022-03-26T18:32:18.534534000', '2022-03-31T18:50:48.087361000',
       '2022-04-02T18:38:26.196010000', '2022-04-04T18:26:53.904150000',
       '2022-04-09T18:44:41.660779000', '2022-04-11T18:32:20.904923000',
       '2022-04-16T18:50:53.974035000', '2022-04-18T18:38:31.913586000',
       '2022-04-20T18:26:57.478599000', '2022-04-25T18:44:39.910946000',
       '2022-04-27T18:32:17.063038000', '2022-05-02T18:50:52.394776000',
       '2022-05-04T18:38:32.304289000', '2022-05-06T18:26:59.822251000',
       '2022-05-11T18:44:47.359594000', '2022-05-13T18:32:26.621212000',
       '2022-05-18T18:51:00.436500000', '2022-05-20T18:38:38.910059000',
       '2022-05-22T18:27:04.890363000', '2022-05-27T18:44:48.250613000',
       '2022-05-29T18:32:28.727749000', '2022-06-03T18:51:05.898078000',
       '2022-06-05T18:38:45.960055000', '2022-06-07T18:27:13.702518000',
       '2022-06-12T18:45:01.955709000', '2022-06-14T18:32:41.535010000',
       '2022-06-19T18:51:16.329680000', '2022-06-21T18:38:55.358354000',
       '2022-06-23T18:27:22.024943000', '2022-06-28T18:45:07.537274000',
       '2022-06-30T18:32:46.074409000', '2022-07-05T18:51:18.547117000',
       '2022-07-07T18:38:56.758064000', '2022-07-09T18:27:22.619818000',
       '2022-07-14T18:45:06.061210000', '2022-07-16T18:32:45.301782000',
       '2022-07-21T18:51:22.298885000', '2022-07-23T18:39:02.212774000',
       '2022-07-25T18:27:29.794387000', '2022-07-30T18:45:17.725644000',
       '2022-08-01T18:32:57.253868000', '2022-08-06T18:51:32.302026000',
       '2022-08-08T18:39:11.550933000', '2022-08-10T18:27:38.463263000',
       '2022-08-15T18:45:24.556392000', '2022-08-17T18:33:03.284411000',
       '2022-08-22T18:51:36.142715000', '2022-08-24T18:39:14.519012000',
       '2022-08-26T18:27:40.622198000', '2022-08-31T18:45:24.812984000',
       '2022-09-02T18:33:02.829160000', '2022-09-07T18:51:38.534099000',
       '2022-09-09T18:39:17.982204000', '2022-09-16T18:45:31.348066000',
       '2022-09-18T18:33:10.105949000', '2022-09-23T18:51:42.960724000',
       '2022-09-25T18:39:21.265191000', '2022-09-27T18:27:47.207724000',
       '2022-10-02T18:45:30.599018000', '2022-10-04T18:33:08.094669000',
       '2022-10-09T18:51:37.077805000', '2022-10-11T18:39:13.505811000',
       '2022-10-13T18:27:37.669230000', '2022-10-18T18:45:26.431834000',
       '2022-10-20T18:33:06.096236000', '2022-10-25T18:51:41.069181000',
       '2022-10-27T18:39:20.161786000', '2022-10-29T18:27:46.833208000',
       '2022-11-03T18:45:31.741431000', '2022-11-05T18:33:09.719987000',
       '2022-11-10T18:51:39.986856000', '2022-11-12T18:39:17.046216000',
       '2022-11-14T18:27:41.667336000', '2022-11-19T18:45:24.533861000',
       '2022-11-21T18:33:03.664615000', '2022-11-26T18:51:37.421591000',
       '2022-11-28T18:39:16.031024000', '2022-11-30T18:27:42.274155000',
       '2022-12-05T18:45:50.377394000', '2022-12-07T18:33:04.379724000',
       '2022-12-12T18:51:34.594981000', '2022-12-14T18:39:11.696450000',
       '2022-12-16T18:27:36.338531000', '2022-12-21T18:45:40.193300000',
       '2022-12-23T18:32:54.280347000', '2022-12-28T18:51:28.515857000',
       '2022-12-30T18:39:31.118057000', '2023-01-01T18:27:33.470203000',
       '2023-01-06T18:46:05.172439000', '2023-01-08T18:32:54.991488000',
       '2023-01-15T18:39:00.610083000', '2023-01-17T18:27:48.368126000',
       '2023-01-22T18:45:13.659001000', '2023-01-24T18:32:53.086260000',
       '2023-01-29T18:51:26.822968000', '2023-01-31T18:39:05.296909000'],
      dtype='datetime64[ns]')
    • y
      (y)
      float64
      1.413e+07 1.413e+07 ... 1.339e+07
      units :
      metre
      resolution :
      -30.0
      crs :
      epsg:32650
      array([14126265., 14126235., 14126205., ..., 13390995., 13390965., 13390935.])
    • x
      (x)
      float64
      5.1e+06 5.1e+06 ... 5.835e+06
      units :
      metre
      resolution :
      30.0
      crs :
      epsg:32650
      array([5099805., 5099835., 5099865., ..., 5834535., 5834565., 5834595.])
    • spatial_ref
      ()
      int32
      32650
      spatial_ref :
      PROJCS["WGS 84 / UTM zone 50N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",117],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32650"]]
      grid_mapping_name :
      transverse_mercator
      array(32650, dtype=int32)
    • qa_pixel
      (time, y, x)
      uint16
      dask.array<chunksize=(2, 4903, 3072), meta=np.ndarray>
      units :
      bit_index
      nodata :
      1
      flags_definition :
      {'snow': {'bits': 5, 'values': {'0': 'not_high_confidence', '1': 'high_confidence'}}, 'clear': {'bits': 6, 'values': {'0': 'not_clear', '1': 'clear'}}, 'cloud': {'bits': 3, 'values': {'0': 'not_high_confidence', '1': 'high_confidence'}}, 'water': {'bits': 7, 'values': {'0': 'land_or_cloud', '1': 'water'}}, 'cirrus': {'bits': 2, 'values': {'0': 'not_high_confidence', '1': 'high_confidence'}}, 'nodata': {'bits': 0, 'values': {'0': False, '1': True}}, 'qa_pixel': {'bits': [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15], 'values': {'1': 'Fill', '2': 'Dilated Cloud', '4': 'Cirrus', '8': 'Cloud', '16': 'Cloud Shadow', '32': 'Snow', '64': 'Clear', '128': 'Water', '256': 'Cloud Confidence low bit', '512': 'Cloud Confidence high bit', '1024': 'Cloud Shadow Confidence low bit', '2048': 'Cloud Shadow Confidence high bit', '4096': 'Snow Ice Confidence low bit', '8192': 'Snow Ice Confidence high bit', '16384': 'Cirrus Confidence low bit', '32768': 'Cirrus Confidence high bit'}, 'description': 'Level 2 pixel quality'}, 'cloud_shadow': {'bits': 4, 'values': {'0': 'not_high_confidence', '1': 'high_confidence'}}, 'dilated_cloud': {'bits': 1, 'values': {'0': 'not_dilated', '1': 'dilated'}}, 'cloud_confidence': {'bits': [8, 9], 'values': {'0': 'none', '1': 'low', '2': 'medium', '3': 'high'}}, 'cirrus_confidence': {'bits': [14, 15], 'values': {'0': 'none', '1': 'low', '2': 'reserved', '3': 'high'}}, 'snow_ice_confidence': {'bits': [12, 13], 'values': {'0': 'none', '1': 'low', '2': 'reserved', '3': 'high'}}, 'cloud_shadow_confidence': {'bits': [10, 11], 'values': {'0': 'none', '1': 'low', '2': 'reserved', '3': 'high'}}}
      crs :
      epsg:32650
      grid_mapping :
      spatial_ref
      Array Chunk
      Bytes 125.25 GiB 57.46 MiB
      Shape (112, 24512, 24494) (2, 4903, 3072)
      Dask graph 2240 chunks in 2 graph layers
      Data type uint16 numpy.ndarray
      		24494 24512 112
    • red
      (time, y, x)
      uint16
      dask.array<chunksize=(2, 4903, 3072), meta=np.ndarray>
      units :
      reflectance
      nodata :
      0
      crs :
      epsg:32650
      grid_mapping :
      spatial_ref
      Array Chunk
      Bytes 125.25 GiB 57.46 MiB
      Shape (112, 24512, 24494) (2, 4903, 3072)
      Dask graph 2240 chunks in 2 graph layers
      Data type uint16 numpy.ndarray
      		24494 24512 112
    • nir08
      (time, y, x)
      uint16
      dask.array<chunksize=(2, 4903, 3072), meta=np.ndarray>
      units :
      reflectance
      nodata :
      0
      crs :
      epsg:32650
      grid_mapping :
      spatial_ref
      Array Chunk
      Bytes 125.25 GiB 57.46 MiB
      Shape (112, 24512, 24494) (2, 4903, 3072)
      Dask graph 2240 chunks in 2 graph layers
      Data type uint16 numpy.ndarray
      		24494 24512 112
    • time
      PandasIndex
      PandasIndex(DatetimeIndex(['2022-02-04 18:45:01.040642', '2022-02-06 18:32:39.001148',
               '2022-02-11 18:51:09.926245', '2022-02-13 18:38:47.450521',
               '2022-02-15 18:27:12.621325', '2022-02-20 18:44:54.824298',
               '2022-02-22 18:32:33.208594', '2022-02-27 18:51:05.324304',
               '2022-03-01 18:38:43.403434', '2022-03-03 18:27:09.162918',
               ...
               '2022-12-30 18:39:31.118057', '2023-01-01 18:27:33.470203',
               '2023-01-06 18:46:05.172439', '2023-01-08 18:32:54.991488',
               '2023-01-15 18:39:00.610083', '2023-01-17 18:27:48.368126',
               '2023-01-22 18:45:13.659001', '2023-01-24 18:32:53.086260',
               '2023-01-29 18:51:26.822968', '2023-01-31 18:39:05.296909'],
              dtype='datetime64[ns]', name='time', length=112, freq=None))
    • y
      PandasIndex
      PandasIndex(Float64Index([14126265.0, 14126235.0, 14126205.0, 14126175.0, 14126145.0,
              14126115.0, 14126085.0, 14126055.0, 14126025.0, 14125995.0,
              ...
              13391205.0, 13391175.0, 13391145.0, 13391115.0, 13391085.0,
              13391055.0, 13391025.0, 13390995.0, 13390965.0, 13390935.0],
             dtype='float64', name='y', length=24512))
    • x
      PandasIndex
      PandasIndex(Float64Index([5099805.0, 5099835.0, 5099865.0, 5099895.0, 5099925.0, 5099955.0,
              5099985.0, 5100015.0, 5100045.0, 5100075.0,
              ...
              5834325.0, 5834355.0, 5834385.0, 5834415.0, 5834445.0, 5834475.0,
              5834505.0, 5834535.0, 5834565.0, 5834595.0],
             dtype='float64', name='x', length=24494))
  • crs :
    epsg:32650
    grid_mapping :
    spatial_ref

Very marginal increase in memory per chunk but the tasks have dropped from 5526 to 4661. Note that this occurs for every measurement and operation so the benefit is significant.

In [20]:
ndvi_unweighted
Out[20]:
<xarray.DataArray (season: 4, y: 24512, x: 24494)>
dask.array<stack, shape=(4, 24512, 24494), dtype=float64, chunksize=(1, 4903, 3072), chunktype=numpy.ndarray>
Coordinates:
  * y            (y) float64 1.413e+07 1.413e+07 ... 1.339e+07 1.339e+07
  * x            (x) float64 5.1e+06 5.1e+06 5.1e+06 ... 5.835e+06 5.835e+06
    spatial_ref  int32 32650
  * season       (season) object 'DJF' 'JJA' 'MAM' 'SON'
xarray.DataArray
  • season: 4
  • y: 24512
  • x: 24494
  • dask.array<chunksize=(1, 4903, 3072), meta=np.ndarray>
    Array Chunk
    Bytes 17.89 GiB 114.91 MiB
    Shape (4, 24512, 24494) (1, 4903, 3072)
    Dask graph 160 chunks in 32 graph layers
    Data type float64 numpy.ndarray
    		24494 24512 4
    • y
      (y)
      float64
      1.413e+07 1.413e+07 ... 1.339e+07
      units :
      metre
      resolution :
      -30.0
      crs :
      epsg:32650
      array([14126265., 14126235., 14126205., ..., 13390995., 13390965., 13390935.])
    • x
      (x)
      float64
      5.1e+06 5.1e+06 ... 5.835e+06
      units :
      metre
      resolution :
      30.0
      crs :
      epsg:32650
      array([5099805., 5099835., 5099865., ..., 5834535., 5834565., 5834595.])
    • spatial_ref
      ()
      int32
      32650
      spatial_ref :
      PROJCS["WGS 84 / UTM zone 50N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",117],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32650"]]
      grid_mapping_name :
      transverse_mercator
      array(32650, dtype=int32)
    • season
      (season)
      object
      'DJF' 'JJA' 'MAM' 'SON'
      array(['DJF', 'JJA', 'MAM', 'SON'], dtype=object)
    • y
      PandasIndex
      PandasIndex(Float64Index([14126265.0, 14126235.0, 14126205.0, 14126175.0, 14126145.0,
              14126115.0, 14126085.0, 14126055.0, 14126025.0, 14125995.0,
              ...
              13391205.0, 13391175.0, 13391145.0, 13391115.0, 13391085.0,
              13391055.0, 13391025.0, 13390995.0, 13390965.0, 13390935.0],
             dtype='float64', name='y', length=24512))
    • x
      PandasIndex
      PandasIndex(Float64Index([5099805.0, 5099835.0, 5099865.0, 5099895.0, 5099925.0, 5099955.0,
              5099985.0, 5100015.0, 5100045.0, 5100075.0,
              ...
              5834325.0, 5834355.0, 5834385.0, 5834415.0, 5834445.0, 5834475.0,
              5834505.0, 5834535.0, 5834565.0, 5834595.0],
             dtype='float64', name='x', length=24494))
    • season
      PandasIndex
      PandasIndex(Index(['DJF', 'JJA', 'MAM', 'SON'], dtype='object', name='season'))

Total task count is sub 100_000 so should be okay but task graph optimisation will take a while. Resulting array is a bit big for the notebook node as stated previously.

The shape spatially is y:15686, x:13707. Standard plots aren't going to work very well for visualising the result in the notebook and the result uses a fair amount of memory so we'll need a different approach.

For now let's visualize the same ROI as the small area before. We stashed that ROI in x_slice, y_slice.

If you haven't already, open the dask dashboard so you can watch the cluster make progress

The code to do this visualisation is basically the same as before except now we specify a slice

In [21]:
%%time
ndvi_unweighted.sel(season='DJF', x=x_slice, y=y_slice).compute().plot(robust=True, size=6, aspect='equal')

CPU times: user 9.6 s, sys: 1.27 s, total: 10.9 s
Wall time: 3min 17s
Out[21]:
<matplotlib.collections.QuadMesh at 0x7f6d03525bd0>

The computation time is relatively short since we are only materialising the result for a subset of the overall dataset.

Visualising all of the data¶

To visualise all of the data we will make use of the dask cluster and some dask-aware visualisation capabilties from holoviews and datashader python libraries. These libraries provide an interactive visualisation capability that leaves the large datasets on the cluster and transmits only the final visualisation to the Jupyter notebook. This is done on the fly so the user can zoom and pan about the dataset in all dimensions and the dask cluster will scale data to fit in the viewport automatically. Details of how this is done and advanced features available is beyond the scope of this dask and ODC course but the manuals are extensive and the basic example here both powerful and useful.

Tip: The datashader pipeline page provides an excellent summary of what's going on.

compute() and persist()¶

The first thing we will do is persist() the results of our calculation to the cluster. This will materialise the results but will keep the result on the cluster (so all lazy tasks are calculated, just like compute() but data locality remains on the cluster). This will ensure the result is readily available for the visualisation. The cluster has plenty of (distributed) memory so there is no reason not to materialise the result on the cluster.

persist() is non-blocking so will return just as soon as task graph optimisation (which is performed in the notebook kernel) is complete. Run the next cell and you will see it takes a few seconds to do task graph optimisation, and once that is complete the Jupyter notebook will be available for use again. At the same time the dask dashboard will show tasks running as the result is computed and left on the cluster.

In [22]:
%%time
on_cluster_result = ndvi_unweighted.persist()
# wait(on_cluster_result)
on_cluster_result

CPU times: user 1.78 s, sys: 44.1 ms, total: 1.83 s
Wall time: 1.8 s
Out[22]:
<xarray.DataArray (season: 4, y: 24512, x: 24494)>
dask.array<stack, shape=(4, 24512, 24494), dtype=float64, chunksize=(1, 4903, 3072), chunktype=numpy.ndarray>
Coordinates:
  * y            (y) float64 1.413e+07 1.413e+07 ... 1.339e+07 1.339e+07
  * x            (x) float64 5.1e+06 5.1e+06 5.1e+06 ... 5.835e+06 5.835e+06
    spatial_ref  int32 32650
  * season       (season) object 'DJF' 'JJA' 'MAM' 'SON'
xarray.DataArray
  • season: 4
  • y: 24512
  • x: 24494
  • dask.array<chunksize=(1, 4903, 3072), meta=np.ndarray>
    Array Chunk
    Bytes 17.89 GiB 114.91 MiB
    Shape (4, 24512, 24494) (1, 4903, 3072)
    Dask graph 160 chunks in 1 graph layer
    Data type float64 numpy.ndarray
    		24494 24512 4
    • y
      (y)
      float64
      1.413e+07 1.413e+07 ... 1.339e+07
      units :
      metre
      resolution :
      -30.0
      crs :
      epsg:32650
      array([14126265., 14126235., 14126205., ..., 13390995., 13390965., 13390935.])
    • x
      (x)
      float64
      5.1e+06 5.1e+06 ... 5.835e+06
      units :
      metre
      resolution :
      30.0
      crs :
      epsg:32650
      array([5099805., 5099835., 5099865., ..., 5834535., 5834565., 5834595.])
    • spatial_ref
      ()
      int32
      32650
      spatial_ref :
      PROJCS["WGS 84 / UTM zone 50N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",117],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32650"]]
      grid_mapping_name :
      transverse_mercator
      array(32650, dtype=int32)
    • season
      (season)
      object
      'DJF' 'JJA' 'MAM' 'SON'
      array(['DJF', 'JJA', 'MAM', 'SON'], dtype=object)
    • y
      PandasIndex
      PandasIndex(Float64Index([14126265.0, 14126235.0, 14126205.0, 14126175.0, 14126145.0,
              14126115.0, 14126085.0, 14126055.0, 14126025.0, 14125995.0,
              ...
              13391205.0, 13391175.0, 13391145.0, 13391115.0, 13391085.0,
              13391055.0, 13391025.0, 13390995.0, 13390965.0, 13390935.0],
             dtype='float64', name='y', length=24512))
    • x
      PandasIndex
      PandasIndex(Float64Index([5099805.0, 5099835.0, 5099865.0, 5099895.0, 5099925.0, 5099955.0,
              5099985.0, 5100015.0, 5100045.0, 5100075.0,
              ...
              5834325.0, 5834355.0, 5834385.0, 5834415.0, 5834445.0, 5834475.0,
              5834505.0, 5834535.0, 5834565.0, 5834595.0],
             dtype='float64', name='x', length=24494))
    • season
      PandasIndex
      PandasIndex(Index(['DJF', 'JJA', 'MAM', 'SON'], dtype='object', name='season'))

The on_cluster_result will continue to show as a dask array on the cluster - not actual results. Think of it as a handle that links the Jupyter client to the result on the dask cluster.

The cluster will start the computation, but we can continue working in the notebook. Let's import some new visualization libraries: holoviews and datashader. We'll be using datashader.rasterize to handle the scaling of the full dataset which has many more pixels than what what is being displayed in the notebook.

Notice also there are no bounds set on the dataset, we are viewing the entire result, including the season dimension. holoviews will provide an interface for pan, zoom and time (season) selection and you can use the mouse to move around the data.

Tip: Keep watching your Dask dashboard to see how the calculations are progressing.

In [23]:
import holoviews as hv
import xarray as xr
from holoviews import opts
from holoviews.operation.datashader import rasterize

hv.extension("bokeh", width=800)

holoviews expects the xarray.DataArray to have a name so let's give it one via the name attribute

In [24]:
on_cluster_result.name = "ndvi"
on_cluster_result
Out[24]:
<xarray.DataArray 'ndvi' (season: 4, y: 24512, x: 24494)>
dask.array<stack, shape=(4, 24512, 24494), dtype=float64, chunksize=(1, 4903, 3072), chunktype=numpy.ndarray>
Coordinates:
  * y            (y) float64 1.413e+07 1.413e+07 ... 1.339e+07 1.339e+07
  * x            (x) float64 5.1e+06 5.1e+06 5.1e+06 ... 5.835e+06 5.835e+06
    spatial_ref  int32 32650
  * season       (season) object 'DJF' 'JJA' 'MAM' 'SON'
xarray.DataArray
'ndvi'
  • season: 4
  • y: 24512
  • x: 24494
  • dask.array<chunksize=(1, 4903, 3072), meta=np.ndarray>
    Array Chunk
    Bytes 17.89 GiB 114.91 MiB
    Shape (4, 24512, 24494) (1, 4903, 3072)
    Dask graph 160 chunks in 1 graph layer
    Data type float64 numpy.ndarray
    		24494 24512 4
    • y
      (y)
      float64
      1.413e+07 1.413e+07 ... 1.339e+07
      units :
      metre
      resolution :
      -30.0
      crs :
      epsg:32650
      array([14126265., 14126235., 14126205., ..., 13390995., 13390965., 13390935.])
    • x
      (x)
      float64
      5.1e+06 5.1e+06 ... 5.835e+06
      units :
      metre
      resolution :
      30.0
      crs :
      epsg:32650
      array([5099805., 5099835., 5099865., ..., 5834535., 5834565., 5834595.])
    • spatial_ref
      ()
      int32
      32650
      spatial_ref :
      PROJCS["WGS 84 / UTM zone 50N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",117],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32650"]]
      grid_mapping_name :
      transverse_mercator
      array(32650, dtype=int32)
    • season
      (season)
      object
      'DJF' 'JJA' 'MAM' 'SON'
      array(['DJF', 'JJA', 'MAM', 'SON'], dtype=object)
    • y
      PandasIndex
      PandasIndex(Float64Index([14126265.0, 14126235.0, 14126205.0, 14126175.0, 14126145.0,
              14126115.0, 14126085.0, 14126055.0, 14126025.0, 14125995.0,
              ...
              13391205.0, 13391175.0, 13391145.0, 13391115.0, 13391085.0,
              13391055.0, 13391025.0, 13390995.0, 13390965.0, 13390935.0],
             dtype='float64', name='y', length=24512))
    • x
      PandasIndex
      PandasIndex(Float64Index([5099805.0, 5099835.0, 5099865.0, 5099895.0, 5099925.0, 5099955.0,
              5099985.0, 5100015.0, 5100045.0, 5100075.0,
              ...
              5834325.0, 5834355.0, 5834385.0, 5834415.0, 5834445.0, 5834475.0,
              5834505.0, 5834535.0, 5834565.0, 5834595.0],
             dtype='float64', name='x', length=24494))
    • season
      PandasIndex
      PandasIndex(Index(['DJF', 'JJA', 'MAM', 'SON'], dtype='object', name='season'))
In [25]:
aspect = on_cluster_result.sizes['y']/on_cluster_result.sizes['x']
width = 700
height = int(width*aspect)

ndvi_seasonal_mean_ds = on_cluster_result.to_dataset(
    name="ndvi_seasonal_mean"
)  # holoviews works better with datasets so let's convert the xarray DataArray holding ndvi into a Dataset

The next cell will display the result - when its ready. The rasterize function will calculate a representative pixel for display from the full array on the dask cluster. If you monitor the dashboard you will see small bursts of activity across the workers and quite some waiting whilst data transfers occur to bring all the summary information back and transmit it to the Jupyter notebook.

You can use the controls on the right to pan and zoom around the full image. If you zoom in, rasterize will take a moment to generate a new summary for the current zoom level and show more or less detail. Similarly for panning.

In [26]:
hv_ds = hv.Dataset(ndvi_seasonal_mean_ds)

rasterize(
    hv_ds.to(hv.Image, ["x", "y"], "ndvi_seasonal_mean", ["season"]).opts(
        title="NDVI Seasonal Mean",
        cmap="RdYlGn", # NDVI more green the larger the value. 
        clim=(-0.3, 0.8), # we'll clamp the range for visualisation to enhance the visualisation
        colorbar=True,
        frame_width=width,
        frame_height=height
    )
, precompute=True)
Out[26]:

Be a good dask user - Clean up the cluster resources¶

Disconnecting your client is good practice, but the cluster will still be up so we need to shut it down as well

In [ ]:
client.close()
In [ ]:
cluster.shutdown()
In [ ]: